// Geometric Tools Library
// https://www.geometrictools.com
// Copyright (c) 2022 David Eberly
// Distributed under the Boost Software License, Version 1.0
// https://www.boost.org/LICENSE_1_0.txt
// File Version: 2022.12.16

#pragma once

// Gaussian quadrature estimates the integral of a function f(x) defined on
// [-1,1] using
//    integral_{-1}^{1} f(t) dt = sum_{i=0}^{n-1} c[i]*f(r[i])
// where r[i] are the roots to the Legendre polynomial p(t) of degree n and
//   c[i] = integral_{-1}^{1} prod_{j=0,j!=i} (t-r[j]/(r[i]-r[j]) dt
// To integrate over [a,b], a transformation to [-1,1] is applied
// internally:  x - ((b-a)*t + (b+a))/2. The Legendre polynomials are
// generated by
//   P[0](x) = 1, P[1](x) = x,
//   P[k](x) = ((2*k-1)*x*P[k-1](x) - (k-1)*P[k-2](x))/k, k >= 2
// Implementing the polynomial generation is simple, and a numerical method is
// required to compute polynomial roots. The challenging task is to develop an
// efficient algorithm for computing the coefficients c[i] for a specified
// degree. The degree must be two or larger.

#include <GTL/Mathematics/RootFinders/RootsPolynomial.h>
#include <GTL/Utility/Exceptions.h>
#include <GTL/Utility/TypeTraits.h>
#include <array>
#include <cstddef>
#include <functional>
#include <vector>

namespace gtl
{
    template <typename T>
    class IntgGaussianQuadrature
    {
    public:
        static void ComputeRootsAndCoefficients(size_t degree, size_t maxBisections,
            size_t precision, std::vector<T>& roots, std::vector<T>& coefficients)
        {
            std::vector<Polynomial1<T>> poly(degree + 1);
            poly[0].SetDegree(1);
            poly[0][0] = C_<T>(1);
            poly[1].SetDegree(2);
            poly[1][0] = C_<T>(0);
            poly[1][1] = C_<T>(1);
            for (size_t n = 2; n <= degree; ++n)
            {
                T mult0 = static_cast<T>(n - 1) / static_cast<T>(n);
                T mult1 = static_cast<T>(2 * n - 1) / static_cast<T>(n);
                poly[n].SetDegree(n + 1);
                poly[n][0] = -mult0 * poly[n - 2][0];
                for (size_t i = 1; i <= n - 2; ++i)
                {
                    poly[n][i] = mult1 * poly[n - 1][i - 1] - mult0 * poly[n - 2][i];
                }
                poly[n][n - 1] = mult1 * poly[n - 1][n - 2];
                poly[n][n] = mult1 * poly[n - 1][n - 1];
            }

            Polynomial1<RootsPolynomial::BSN> finalPoly(poly[degree].GetDegree());
            for (size_t i = 0; i <= finalPoly.GetDegree(); ++i)
            {
                Convert(poly[degree][i], finalPoly[i], precision);
            }

            RootsPolynomial finder(maxBisections, precision);
            std::vector<RootsPolynomial::BSN> ratRoots{};
            finder(finalPoly, ratRoots);
            roots.resize(ratRoots.size());
            for (size_t i = 0; i < ratRoots.size(); ++i)
            {
                roots[i] = ratRoots[i];
            }

            coefficients.resize(roots.size());
            size_t n = roots.size() - 1;
            std::vector<T> subroots(n);
            for (size_t i = 0; i < roots.size(); ++i)
            {
                T denominator = C_<T>(1);
                for (size_t j = 0, k = 0; j < roots.size(); ++j)
                {
                    if (j != i)
                    {
                        subroots[k++] = roots[j];
                        denominator *= roots[i] - roots[j];
                    }
                }

                std::array<T, 2> delta =
                {
                    -C_<T>(1) - subroots.back(),
                    +C_<T>(1) - subroots.back()
                };

                std::vector<std::array<T, 2>> weights(n);
                weights[0][0] = C_<T>(1, 2) * delta[0] * delta[0];
                weights[0][1] = C_<T>(1, 2) * delta[1] * delta[1];
                for (size_t k = 1; k < n; ++k)
                {
                    T dk = static_cast<T>(k);
                    T mult = -dk / (dk + C_<T>(2));
                    weights[k][0] = mult * delta[0] * weights[k - 1][0];
                    weights[k][1] = mult * delta[1] * weights[k - 1][1];
                }

                struct Info
                {
                    Info()
                        :
                        numBits(0),
                        product{ C_<T>(0), C_<T>(0) }
                    {
                    }

                    size_t numBits;
                    std::array<T, 2> product;
                };

                size_t numElements = (static_cast<size_t>(1) << (n - 1));
                std::vector<Info> info(numElements);
                info[0].numBits = 0;
                info[0].product[0] = C_<T>(1);
                info[0].product[1] = C_<T>(1);
                for (size_t ipow = 1, r = 0; ipow < numElements; ipow <<= 1, ++r)
                {
                    info[ipow].numBits = 1;
                    info[ipow].product[0] = -C_<T>(1) - subroots[r];
                    info[ipow].product[1] = +C_<T>(1) - subroots[r];
                    for (size_t m = 1, j = ipow + 1; m < ipow; ++m, ++j)
                    {
                        info[j].numBits = info[m].numBits + 1;
                        info[j].product[0] =
                            info[ipow].product[0] * info[m].product[0];
                        info[j].product[1] =
                            info[ipow].product[1] * info[m].product[1];
                    }
                }

                std::vector<std::array<T, 2>> sum(n);
                std::array<T, 2> zero2 = { C_<T>(0), C_<T>(0) };
                std::fill(sum.begin(), sum.end(), zero2);
                for (size_t k = 0; k < info.size(); ++k)
                {
                    sum[info[k].numBits][0] += info[k].product[0];
                    sum[info[k].numBits][1] += info[k].product[1];
                }

                std::array<T, 2> total = zero2;
                for (size_t k = 0; k < n; ++k)
                {
                    total[0] += weights[n - 1 - k][0] * sum[k][0];
                    total[1] += weights[n - 1 - k][1] * sum[k][1];
                }

                coefficients[i] = (total[1] - total[0]) / denominator;
            }
        }

        // Gaussian quadrature integration.
        static T Integrate(std::vector<T> const& roots,
            std::vector<T>const& coefficients, T const& a, T const& b,
            std::function<T(T)> const& integrand)
        {
            T radius = C_<T>(1, 2) * (b - a);
            T center = C_<T>(1, 2) * (b + a);
            T result = C_<T>(0);
            for (size_t i = 0; i < roots.size(); ++i)
            {
                result += coefficients[i] * integrand(radius * roots[i] + center);
            }
            result *= radius;
            return result;
        }

    private:
        template <typename SourceType = T, IsFPType<SourceType> = 0>
        inline static void Convert(SourceType const& source,
            RootsPolynomial::BSN& target, size_t = 0)
        {
            target = RootsPolynomial::BSN(source);
        }

        template <typename SourceType = T, IsAPType<SourceType> = 0>
        inline static void Convert(SourceType const& source,
            RootsPolynomial::BSN& target, size_t precision)
        {
            gtl::Convert(source, precision, APRoundingMode::TO_NEAREST, target);
        }

    private:
        friend class UnitTestIntgGaussianQuadrature;
    };
}
